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We discuss the application of the local lattice technique of Maggs and Rossetto |l| to 
problems that involve the motion of objects with different dielectric constants than the 
background. In these systems the simulation method produces a spurious interaction force 
which causes the particles to move in an unphysical manner. We show that this term can 
be removed using a variant of a method known from high-energy physics simulations, the 
multiboson method, and demonstrate the effectiveness of this corrective method on a system 
of neutral particles. We then apply our method to a one-component plasma to show the 
effect of the spurious interaction term on a charged system. 

I. INTRODUCTION 

The long-range character of the electrostatic Coulomb interaction lies at the root of the com- 
putational difficulties encountered in the simulation of many systems of biophysical interest, in 
which one wishes to understand the thermodynamics of a Coulomb gas of charged particles - 
or macroions — moving in an environment of spatially varying dielectric. Recently, Maggs and 
collaborators have suggested rewriting the problem of a classical Coulomb gas in a local 

lattice framework in which each charged particle responds only to the local electric field, which is 
also updated so that Gauss' Law is respected at each point in the simulation. Most applications 
of the original local algorithm, as well as a variety of suggested improvements P, have dealt 
primarily with situations in which the dielectric constant is not spatially varying. In this case, 
the transverse (or "curl") part of the electric field, which is unconstrained in the partition sum, 
decouples from the physics of the charged ions, so that the method correctly calculates the classical 
electrostatic energy of the charges. 

In this paper, we consider situations in which the dielectric constant becomes dynamical: i.e. is 
spatially varying, and in a way that depends on the location of the charges in the system, so that 
it changes in the course of the Monte Carlo simulation. In a molecular dynamics simulation, for 
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example, one effectively determines the electrostatic energy of the system at each step by solving 
the appropriate Poisson equation (in the presence of the varying dielectric) . Of course, the resulting 
energy expression still suffers from the problem of having to include contributions from all pairs of 
charged particles in the system. In the local lattice approach, on the other hand, the unconstrained 
transverse part of the functional integral over the electric field no longer decouples from the charged 
particle dynamics, and leads |2| to an apparent attractive dipole-dipole interaction, even when the 
system is completely electrically neutral, with zero classical electrostatic energy. Although such 
interactions certainly exist at a quantum level (e.g. Casimir, Keesom, van der Waals forces), our 
interest here is in a purely classical calculation. 

Our purpose in this paper is to show that the local algorithm of Maggs et al can be extended 
by addition of a set of boson fields — with a completely local hamiltonian — which remove any 
unphysical contributions from the curl part of the electric field, so that the thermodynamics reflects 
the classical partition function of a charged system with electrostatic energy H es = J as 
desired. In Section[n]we describe in detail the structure of the transverse field contribution in a local 
lattice formulation of the Coulomb gas. In particular, we derive the form of the determinant induced 
by this contribution. In Section ITlTl we describe the use of a multiboson technique familiar in lattice 
gauge theory ( "unquenched" quantum chromo dynamics) to eliminate the spurious determinant 
factor. In Section llVl we test the method by calculating the density-density structure factor, both 
with a system of neutral dielectric particles and with a one-component charged plasma. Finally, 
in Section we summarize our results and discuss the outlook for further applications. 



II. TRANSVERSE MODE CONTRIBUTIONS IN SYSTEMS WITH VARYING 

DIELECTRIC 

In this section we shall show how the local algorithm of references can be 

modified to prevent the appearance of spurious dipole-like interactions that are not present in the 
classical electrostatics of the system when the dielectric constant varies both spatially and in the 
course of the simulation (e.g. if there are mobile entities with dielectric different from that of the 
ambient medium). We illustrate the point first in the case of a continuous system. Later, a spatial 
lattice will be introduced to make the functional integrations precise. 

Consider the partition function of a system consisting of N free charges (mobile or fixed) e% at 
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locations r*j, thereby producing a free charge density 

p{r) =^ei5{r-ri) (1) 

i 

where the system is also described by a linear dielectric function (here assumed isotropic) e(r). 
Note that in general the dielectric function e(r) may depend on the locations of the free charges fi, 
which may be embedded in regions of varying dielectric. This dependence should be kept in mind, 
although for notational convenience it will be suppressed in the following. The electric displacement 
can be broken into longitudinal and transverse parts using the general Helmholtz decomposition, 

D{r) = -e(r)V(p(r) + V x A(r) (2) 
= D\f) + J D tr (r). (3) 

As the transverse and longitudinal components are orthogonal to each other, we have 

dr^=[dr^ + [dr^f (4) 
e(r) J dr e(r) + J e(r) ' [ > 

If the constraint of Maxwell's second law is explicitly imposed the transverse part of the electric 

displacement vanishes and the electrostatic energy of the system is given by 

1 r DW(r) 2 

and the canonical partition function for the system at inverse temperature becomes 

r N 

Z = / II dr ie -P H ™ (6) 

i=l 

where i5" must be determined by first solving —V • (eV0) = p, from which one obtains Z)" = 
— e(r)V0(r). Note that the (nonexistent) transverse part of the displacement field plays no role in 
this result. 

On the other hand, the partition function proposed in P] includes an integral over both the 
transverse and longitudinal parts of the electric displacement and reads simply 

N 



Z' = JfldfiYlVDi^diV-D-pi^e-^I^ 32 /^ (7) 

i=l r 

r N 

= / ]Jdfi Y[VD^(f)VD tv (r)5(V ■ DW - p(r)) 

i=l r 

Xe"! / dr&\(r) 2 /e(r) e -§ J drD^(r) 2 /e(r) ^ 

as a result of the identity Eq. 0J It is apparent that the integration over transverse degrees of 
freedom D tr in Eq. |H] necessarily introduces a dependence on the particle locations r, through the 
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dependence on e(r) absent in the conventional electrostatic energy of Eq. |SJ The exact form of this 
spurious e-dependence can be uncovered by considering the form of Z' for the special case ej = of 
uncharged particles, where D" = and the entire dependence on e(r) derives from the transverse 
part of the functional integral in Z' . 

Of course, in the real world quantum fluctuations of the transverse (and longitudinal) parts of 
the field do exist, and would have to be included in a fully quantum-mechanical treatment of the 
Coulomb gas problem. This is not, however, the problem being addressed here, where we are com- 
puting the purely classical partition function of a classical Coulomb gas. Accordingly, the spurious 
interaction potentials induced by the integration over transverse degrees of freedom in Eq. [HI should 
not be identified with Keesom potentials, for example, which have a quantum mechanical origin, 
and a strength dependent on Planck's constant, which appears nowhere in our classical partition 
function. Such non-Coulombic potentials can of course be included phenomenologically in our 
classical treatment as an explicit separate contribution to the energy function (with an appropri- 
ate interaction strength). For example, hydrophobic interactions, which depend on intermolecular 
forces and thus are quantum mechanical in origin, play a significant role in the organization and 
function of molecular level biophysical structures (e.g., membranes, proteins and nucleic acids). 
Hydrophobic effects can be included empirically by adding to the electrostatic free energy a term 
proportional (with constant of proportionality often called the "surface tension") to the surface 
area of the membrane, protein, or nucleic acid that is exposed to water [3, S|- 

Before investigating the special case of neutral particles, we shall go over to a lattice discretiza- 
tion of the system. We imagine a lattice Coulomb gas with displacement vector field D nfJi defined 
on lattice links, where lattice sites are denoted n and \x = 1,2,3 indicating the spatial direction 
of the link. Likewise, it turns out to be convenient to associate a dielectric value with each link 
(rather than site), so that the dielectric function on the lattice b ecomes For example, in prob- 
lems involving macroions extending over several lattice sites and with dielectric constant differing 
from that of the ambient medium, links crossing the surface of the macroion can be chosen to have 
an appropriately interpolated value of the dielectric constant. Alternatively, this can be regarded 
as a generalization to allow nonisotropic systems where the principal axes of the dielectric tensor 
coincide with the lattice axes. The entire discussion given below can readily be generalized to the 
case of a completely general dielectric tensor field. Introducing left (resp. right) lattice derivatives 
A (resp. A), the Helmholtz decomposition of the displacement field on the lattice becomes 

-D n /i — ^n/i^fM^n 5 £fj,va Ajy-Ajj^ (9) 
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where the electrostatic potential (lattice site field) (ft n satisfies the Poisson equation 

- Y A M (e nM A M n ) = p n (10) 

i 

In the absence of free charges, = 0, the local form of the partition function Z' becomes 

N 

Z'(e t = 0) = J J] dr^(e) (12) 

i=l 

where the e-dependence is entirely due to the transverse degrees of freedom and enters through the 
function 

Hz) = /n^ 5 (E VV)e~^ E ^ D ^ AnA1 (13) 

nfi fi 

= I "jJdAnn^e^n^^-lEn^A"*. (14) 

n nfi 

= I II dX n II di V e "' E "" ^ D - A ^ (15) 

= CW^rJ Y[dX n e-^^ e ^ {A ^ )2 (16) 
= C'nV^det-^-^A^A^) (17) 

nfi fi 

In going from l|13j) to 1)141) we have introduced an auxiliary field X n to implement the Gauss' 
Law constraint (for an everywhere neutral system), allowing the Gaussian integration over the 
displacement field D nfl to be carried out explicitly. The constants C and C are independent of e 
and of no further interest. The integration over the auxiliary A n field (also Gaussian!) then yields 
the determinant of the indicated operator, whose action on a lattice site field takes the explicit 
form 

6 6 

MX„, = (—^2 A M e At A M )A n = ^ tniK ~ Y € m^n+i (18) 

ft i=l i=l 

where the index i now runs over both positive (i = 1,2,3) and negative (i = 4,5,6) spatial 
directions, and e n j = e n+ i t i-3 for i = 4, 5, 6. 

In order to generate an ensemble of configurations based on the partition function Eq. |§] arising 
from the purely physical electrostatic energy, the spurious functional dependence of J-{e) must be 
removed: in other words, we should insert a factor 

r-\e) = e"^ ^ log ^det+i (M) (19) 
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into the partition integral (J7J) to remove the effect of the transverse field modes. In the event that 
the dielectric function is truly spatially- independent, or is spatially varying but frozen throughout 
the simulation, adding the factor of J r_1 (e) is unnecessary, as the transverse integration decouples 
from the dynamics of the problem. 

The problem that we are faced with here is well known: positive (fractional or integral) powers of 
determinants of local operators are intrinsically nonlocal, unlike negative powers, which may be re- 
expressed as integrals over auxiliary fields with local actions (cf. Eqs. ((T5|) through (|17|)1. In lattice 
quantum chromodynamics (QCD), for example, the inclusion of virtual quark-antiquark processes 
leads to precisely the positive power of the determinant of the Dirac operator in the path integral 
for the system, greatly increasing the difficulty of simulations in the full ( "unquenched" ) version 
of the theory in comparison to the truncated ("quenched") version where the quark determinant 
is simply dropped [13]. 

In the next section we shall show how the multiboson technique introduced by Liischer (|9|) for 
approximately computing the positive power of quark determinants can be used to write a purely 
local Hamiltonian in terms of a set of auxiliary scalars ("multibosons" in the lattice QCD language) 
which is readily susceptible to Monte Carlo simulation and will allow us to implement the correct 
electrostatic partition function (jSJ) for systems with general dielectric makeup. 



III. ELIMINATION OF TRANSVERSE CONTRIBUTIONS USING MULTIBOSON 

FIELDS 

The representation of determinants of local operators by integrals over multiple auxiliary scalar 
fields was first introduced by Liischer in the context of unquenched lattice QCD. The essential 
point is to find a uniform polynominal approximation to the function 1 / s in the interval [5, 11 for 
small 5. In terms of the complex roots of the polynominal = + iv^, a convenient choice 3] is 
the Chebyshev polynominal of order 2Nb with 

1 N B 



P(s) = C H((s - Mfc ) 2 + v\) (20) 
s k=i 

1 ?Trk 

m _-(! + „(! _.»__) (21) 

rz . 2nk 

Vk = V6sm m^TT (22) 
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The representation of Eq. [20] can be extended to the determinant of a real positive symmetric 
operator with spectrum in the interval [0, 1] 

N B 

det+5 (M) ~ J] det ~^ d M ~ ^? + v k) ( 23 ) 
fc=i 

where the representation becomes exact in the limit Nb — * oo, 5 — > 0. However, we shall see that 
accurate results can be obtained with surprisingly small values of iVg (=# of auxiliary scalar fields, 
see below). 

The essential idea of the multiboson approach is to replace the determinant factors on the right 
of Eq. |221by integrals over a set of auxiliary fields 4>n , k = 1,2, ..Nb, with local actions, as follows 




where, once again, C is an irrelevant constant that can henceforth be neglected. As pointed out 
previously, our polynominal representation assumes that the spectrum of the operator M in Eq. 
El lies in the interval [0, 1]. Let us assume that the dielectric function e of our system is bounded 
above by the value eo (frequently, in biophysical simulations, this is ~80, corresponding to the 
dielectric constant of an aqueous ambient medium). Recalling that the lattice Laplacian operator 
has largest eigenvalue equal to 12, one easily shows that the operator 

M = J-(-A^A M ) (25) 

has a spectrum contained in the unit interval provided K > 12. In the multiboson approach, 
the auxiliary scalar field <j>^ is entrusted with resolving the spectrum of the operator A4 in the 
neighborhood of in a region of width i/j. (see |9|). For k near 1 or Nb, is small and /i^ 
is close to zero or one, and only a small region of the spectrum is accurately treated. We can 
eliminate the lack of resolution for finite Nb at the upper end of the spectrum by choosing K 
somewhat larger than 12 (in the simulations reported below, we typically take -fT=13), but if the 
spectrum of M is very dense near the origin, we will necessarily be forced to use a large value of 
Nb- Fortunately, in the systems we have so far simulated, this does not appear to be the case. This 
is a fortunate distinction from the case of lattice quantum chromodynamics, where chiral symmetry 
breaking necessarily implies a dense spectrum of eigenvalues of the quark Dirac operator at the 
origin! Instead, in the systems of concern to us, the region with large dielectric eo occupies almost 
all of the volume of the system, and the spectral density of Ai does not differ appreciably from 
that of the free Laplacian. 
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To summarize our proposal, we shall consider a system of mobile charged entities (charges ej, 
locations fi) with total Hamiltonian 

H tot = H es + V nc (r i ) (26) 

where the purely electrostatic energy H es is given in Eq. while all other non-Coulombic energy 
effects (exclusion, Keesom, van der Waals, hard-core, soft-core, etc. etc.) are included in V nc , 
with appropriate phenomenological values. The latter are not typically problematic as they are 
essentially short-range effects. From Eqs. (|17I19I24(I we find that the correct expression for the 
partition function is 



Z = f J] dfi J] dD nfl 5(Y, K D n» - Pn) 

x^- 1 (e)e- (3Vnc e~^^ D "'" /en ' 1 (27) 



n dn n dD nil n #w $e - Pn ) 

i n^i kn A* 

xe -^-4 E„ P lo s (*»#») e -f E„ M 

• A? 



xe 



E^ fc) ((A4-^) 2 +^ fc) (28) 



where A'i is scaled as in Eq. [23 and the integrals over fi are taken to be sums when the charges 
are constrained to be on lattice sites. The effective Hamiltonian appearing in the exponent terms 
of (|28|) is completely local: to simulate the system, we need only update, in some conveniently 
chosen order, (i) the locations fj of the charged entities, (ii) the lattice displacement field D nfl 
(respecting locally the Gauss' Law constraint), and (iii) the Nb auxiliary scalar fields (fi^ k \ in all 
cases according to the indicated Boltzmann weight. 



IV. RESULTS 



In order to see the effect of the additional dipole-like interaction present in the uncorrected 
simulation and study the effectiveness of our correction scheme, we have simulated a system of 
neutral particles similar to the system studied in Ref. Q|. In classical electrodynamics, neutral 
particles, even those with a different dielectric constant from the background, do not interact and 
therefore the density-density structure factor, S(q), is constant as a function of q. In contrast, the 
simulations in Ref. [2| show that using the local lattice approach (not corrected to remove the extra 
dipole term) to simulate neutral particles with different dielectric constants from the background 
results in a clearly non-constant density-density structure factor. In this section we show that the 
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FIG. 1: Residual from polynominal approximation used in the multiboson correction factor for the simula- 
tions in this paper. 

multiboson method presented in this paper is able to effectively remove the dipole terms so that 
the neutral particle simulation gives a constant density-density structure factor. We then apply the 
method to the more interesting case of a one-component plasma with charged particles of different 
dielectric constant than the background. 

In the simulations described below we have compared the effect of using varying number of 
multiboson fields to determine the sensitivity of our results to the low eigenvalues of the operator 
M discussed in the preceding section. In particular, we have simulated neutral systems with 
-/V#=4 boson fields (and 5 =0.07) and charged systems with Nb =4,6,8 (and 5 =0.07,0.07, and 
0.05 respectively). It is instructive to see the accuracy of the polynominal approximation P(s) to 
1/s for these choices by plotting the residual P(s) — 1/s, as shown in Fig. ^ 

We first consider a system of 1000 neutral particles on a 16 3 lattice. The particles are constrained 
to the lattice site and only one particle is allowed on each lattice site. The background dielectric 
constant was set to 1.0, and the dielectric constant of the particles was set to 0.2, except in 
the uniform dielectric simulations where the particles had dielectric constant 1.0. The dielectric 
constant along a link in the /U direction from a lattice site n is defined through the relation 



where e n and e n +^ are either the background dielectric constant or the particle dielectric constant 
depending on whether there is a particle on site n or site n + /i respectively. The dimensionless 
inverse temperature, j3 = 4-7re 2 /a, was set to 0.25. We performed 5000 Monte Carlo warmup sweeps, 
followed by 40, 000 Monte Carlo sweeps. To update the electric field we used the local heat bath 
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FIG. 2: Density-density structure factor for the neutral system. Shown is the results from both a simulation 
with the multiboson correction factor and an uncorrected simulation. The dotted line gives the desired q 
independant structure factor. 



method discussed in Ref. [5j. To update the neutral particles at each Monte Carlo step we chose 
1000 particles, which we attempt to move using the methods discussed in Section UTTI 

Figure shows the density-density structure factor for both the corrected and uncorrected 
simulation. In the uniform case the structure factor is independent of q and is given by N(l— ), 
where N is the number of particles and V is the volume of the system (the number of points in 
the lattice). The spurious dipole-like interactions present in the uncorrected simulation give a 
density-density structure factor that differs significantly from the desired flat structure factor of 
the uniform case, while the simulation with the multiboson correction reproduces the analytically 
calculated flat structure factor of the neutral system. The acceptance rate for particle moves drops 
from 30% in the uncorrected simulation to 15% in the multiboson simulation. The uncorrected 
simulation took 11 hours to run on a single processor workstation and the corrected simulation 
took 55 hours on a similar workstation. 

After verifying that the multiboson method was able to correctly remove the non-physical dipole 
interaction from the simulation of neutral particles, we applied the multiboson technique to the 
simulation of a one-component charged plasma. The particles in the system are positively charged 
with one electron charge each and the background sites are uniformly given enough negative charge 
so that the entire system is charge neutral. The parameters of the system are the same as in the 
neutral particle case except the lattice is a 32 3 lattice, there are 8000 particles, 8000 particles are 
updated each Monte Carlo step, and the dielectric constant of the particles in the non-uniform case 
was taken to be 0.05. In this case we have also repeated the simulations with varying numbers of 
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qL/2n 

FIG. 3: Density-density structure factor for the charged particles in the one-component plasma. Shown are 
the results from the uncorrected simulation, the simulation with the multiboson correction factor, and a 
simulation where the dielectric constant is uniform. The errorbars are smaller then the symbol size. 

boson fields (N B =A,6,8). 

In Fig. |31 we show the results for the density-density structure function for three cases: the 
fully multiboson corrected structure function for this system, the uncorrected structure function 
switching off multiboson contributions, and for a charged plasma with uniform dielectric (e par t = 
e bg = 1.0). It is apparent from the figure that the removal of the spurious interactions induced by 
the transverse part of the field in the nonuniform dielectric case results in a qualitative modification 
of the shape of the structure function. Also, the results differ significantly between the cases of 
uniform and nonuniform dielectric. This difference may play an important role in the behavior of 
systems from biological and chemical physics, so it is important to be able to reliably calculate the 
effect of having a non-uniform dielectric which changes dynamically in the course of the simulation. 
The results in Fig. |3] were obtained using A r s=4 multiboson fields to estimate the transverse electric 
field determinant. 

Of course, the polynominal approximation Eq . 1201 cannot accurately represent very small eigen- 
values of A4, so it is important to check the sensitivity of our results to the number of boson 
fields used. This comparison, again for the case of the charged one-component plasma on a 32 3 
lattice described above, is shown in Fig. 4. Certainly the results for A r s=4,6, and 8 are in complete 
qualitative agreement over the entire range of q. For larger lattices, one will need to increase iVg 
to deal with the larger dynamical range in the eigenspectrum of M, and, as we saw previously, 
this in turn results in a drop in the acceptance rate for particle moves. A similar problem in the 
use of multiboson fields in unquenched QCD has been addressed by use of a hybrid scheme in 
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FIG. 4: Density-density structure factor for the charged particles in the one-component plasma. Shown are 
the results from simulations corrected with differing numbers of bosonic field in the multiboson correction 
factor. The errorbars are smaller then the symbol size. Only a portion of the larger q results are shown to 
improve the legibility of the plot. 

which the number of multiboson fields is held fixed, and the low eigenvalues of M. treated exactly 
by a Lanczos algorithm. 



Systems from biological and chemical physics frequently have mobile charged elements with 
different dielectric constant from the background. The local lattice approach provides an efficient 
method to simulation these systems. Unfortunately, a spurious force remains as an artifact of the 
simulation method. This force can be efficiently removed using a set of bosonic fields to approximate 
a non-local counteracting force. This method is able to effectively remove the spurious term in 
both neutral and charged systems. 



V. 



CONCLUSION 
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